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Abstract 

This paper presents a geometric- variational approach to continuous and dis- 
crete mechanics and field theories. Using multisymplectic geometry, we show 
that the existence of the fundamental geometric structures as well as their 
preservation along solutions can be obtained directly from the variational prin- 
ciple. In particular, we prove that a unique multisymplectic structure is ob- 
tained by taking the derivative of an action function, and use this structure to 
prove covariant generalizations of conservation of symplecticity and Noether's 
theorem. Natural discretization schemes for PDEs, which have these important 
preservation properties, then follow by choosing a discrete action functional. 
In the case of mechanics, we recover the variational symplectic integrators of 
Veselov type, while for PDEs we obtain covariant spacetime integrators which 
conserve the corresponding discrete multisymplectic form as well as the discrete 
momentum mappings corresponding to symmetries. We show that the usual 
notion of symplecticity along an infinite-dimensional space of fields can be natu- 
rally obtained by making a spacetime split. All of the aspects of our method are 
demonstrated with a nonlinear sine-Gordon equation, including computational 
results and a comparison with other discretization schemes. 
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1 Introduction 

The purpose of this paper is to develop the geometric foundations for multisym- 
plectic- momentum integrators for variational partial differential equations (PDEs). 
These integrators are the PDE generalizations of symplectic integrators that are 
popular for Hamiltonian ODEs (see, for example, the articles in Marsden, Patrick 
and Shadwick [1996], and especially the review article of McLachlan and Scovel 
[1996]) in that they are covariant spacetime integrators which preserve the geometric 
structures of the system. 

Because of the covariance of our method which we shall describe below, the 
resulting integrators are spacetime localizable in the context of hyperbolic PDEs, 
and generalize the notion of symplecticity and symmetry preservation in the context 
of elliptic problems. Herein, we shall primarily focus on spacetime integrators; 
however, we shall remark on the connection of our method with the finite element 
method for elliptic problems, as well as the Gregory and Lin [1991] method in 
optimal control. 

Historically, in the setting of ODEs, there have been many approaches devised for 
constructing symplectic integrators, beginning with the original derivations based on 
generating functions (see de Vogelaere [1956]) and proceeding to symplectic Runge- 
Kutta algorithms, the shake algorithm, and many others. In fact, in many areas 
of molecular dynamics, symplectic integrators such as the Verlet algorithm and 
variants thereof are quite popular, as are symplectic integrators for the integration 
of the solar system. In these domains, integrators that are either symplectic or 
which are adaptations of symplectic integrators, are amongst the most widely used. 

A fundamentally new approach to symplectic integration is that of Veselov [1988], 
[1991] who developed a discrete mechanics based on a discretization of Hamilton's 
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principle. Their method leads in a natural way to symplectic-momentum integrators 
which include the shake and Verlet integrators as special cases (see Wendlandt and 
Marsden [1997]). In addition, Veselov integrators often have amazing properties 
with regard to preservation of integrable structures, as has been shown by Moser 
and Veselov [1991]. This aspect has yet to be exploited numerically, but it seems to 
be quite important. 

The approach we take in this paper is to develop a Veselov-type discretization for 
PDE's in variational form. The relevant geometry for this situation is multisymplec- 
tic geometry (see Gotay, Isenberg, and Marsden [1997] and Marsden and Shkoller 
[1997]) and we develop it in a variational framework. As we have mentioned, this 
naturally leads to multisymplectic- momentum integrators. It is well-known that 
such integrators cannot in general preserve the Hamiltonian exactly (Ge and Mars- 
den [1988]). However, these integrators have, under appropriate circumstances, very 
good energy performance in the sense of the conservation of a nearby Hamiltonian 
up to exponentially small errors, assuming small time steps, due to a result of Neish- 
tadt [1984]. See also Dragt and Finn [1979], and Simo and Gonzales [1993]. This 
is related to backward error analysis; see Sanz-Serna and Calvo [1994], Calvo and 
Hairer [1995], and the recent work of Hyman, Newman and coworkers and references 
therein. It would be quite interesting to develop the links with Neishtadt's analysis 
more thoroughly. 

An important part of our approach is to understand how the symplectic nature 
of the integrators is implied by the variational structure. In this way we are able to 
identify the symplectic and momentum conserving properties after discretizing the 
variational principle itself. Inspired by a paper of Wald [1993], we obtain a formal 
method for locating the symplectic or multisymplectic structures directly from the 
action function and its derivatives. We present the method in the context of ordinary 
Lagrangian mechanics, and apply it to discrete Lagrangian mechanics, and both 
continuous and discrete multisymplectic field theory. While in these contexts our 
variational method merely uncovers the well-known differential-geometric structures, 
our method forms an excellent pedagogical approach to those theories. 

Outline of paper. 

§2 In this section we sketch the three main aspects of our variational approach in 
the familiar context of particle mechanics. We show that the usual symplectic 
2-form on the tangent bundle of the configuration manifold arises naturally 
as the boundary term in the first variational principle. We then show that 
application of d 2 = to the variational principle restricted to the space of 
solutions of the Euler-Lagrange equations produces the familiar concept of 
conservation of the symplectic form; this statement is obtained variationally 
in a non-dynamic context; that is, we do not require an evolutionary flow. 
We then show that if the action function is left invariant by a symmetry 
group, then Noether's theorem follows directly and simply from the variational 
principle as well. 

§3 Here we use our variational approach to construct discretization schemes for 
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mechanics which preserve the discrete symplectic form and the associated dis- 
crete momentum mappings. 

§4 This section defines the three aspects of our variational approach in the multi- 
symplectic field-theoretic setting. Unlike the traditional approach of defining 
the canonical multisymplectic form on the dual of the first jet bundle and then 
pulling back to the Lagrangian side using the covariant Legendre transform, 
we obtain the geometric structure by staying entirely on the Lagrangian side. 
We prove the covariant analogue of the fact that the flow of conservative sys- 
tems consists of symplectic maps; we call this result the multisymplectic form 
formula. After variationally proving a covariant version of Noether's theorem, 
we show that one can use the multisymplectic form formula to recover the 
usual notion of symplecticity of the flow in an infinite-dimensional space of 
fields by making a spacetime split. We demonstrate this machinery using a 
nonlinear wave equation as an example. 

§5 In this section we develop discrete field theories from which the covariant in- 
tegrators follow. We define discrete analogues of the first jet bundle of the 
configuration bundle whose sections are the fields of interest, and proceed to 
define the discrete action sum. We then apply our variational algorithm to this 
discrete action function to produce the discrete Euler-Lagrange equations and 
the discrete multisymplectic forms. As a consequence of our methodology, we 
show that the solutions of the discrete Euler-Lagrange equations satisfy the 
discrete version of the multisymplectic form formula as well as the discrete ver- 
sion of our generalized Noether's theorem. Using our nonlinear wave equation 
example, we develop various multisymplectic-momentum integrators for the 
sine-Gordon equations, and compare our resulting numerical scheme with the 
energy-conserving methods of Li and Vu-Quoc [1995] and Guo, Pascual, Ro- 
driguez, and Vazquez [1986]. Results are presented for long-time simulations 
of kink-antikink solutions for over 5000 soliton collisions. 

§6 This section contains some important remarks concerning the variational inte- 
grator methodology. For example, we discuss integrators for reduced systems, 
the role of grid uniformity, and the interesting connections with the finite- 
element methods for elliptic problems. We also make some comments on 
future work. 

2 Lagrangian Mechanics 

Hamilton's Principle. We begin by recalling a problem going back to Euler, 
Lagrange and Hamilton in the period 1740-1830. Consider an n-dimensional con- 
figuration manifold Q with its tangent bundle TQ. We denote coordinates on Q by 
q % and those on TQ by (q l , q 1 ). Consider a Lagrangian L : TQ — > R. Construct the 
corresponding action functional S on C 2 curves q(t) in Q by integration of L along 
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the tangent to the curve. In coordinate notation, this reads 

S{q(t))= j\[q i {t), d ^{t) \ dt. 



(2.1) 



The action functional depends on a and b, but this is not explicit in the notation. 
Hamilton's principle seeks the curves q(t) for which the functional S is stationary 
under variations of q(t) with fixed endpoints; namely, we seek curves q{t) which 
satisfy 



dS(q(t))-5q(t) 



d_ 



e=0 



S{q e (t)) = 



(2.2) 



for all 5q(t) with Sq(a) = Sq(b) = 0, where q e is a smooth family of curves with 
qo = q and (d/de)\ e =oq € = 5q. Using integration by parts, the calculation for this is 
simply 



dS{q(t))-5q(t) = - 



IX 




dt 


(dL 


d dL\ , 

dtWJ diA 


dL 


yoq 1 ~ 


dq 1 



(2.3) 



The last term in ( ^.3| ) vanishes since Sq(a) = Sq(b) = 0, so that the requirement 
(|2.2|) for S to be stationary yields the Euler- Lagrange equations 



dL d dL _ Q 
dq 1 dt dq 1 



(2.4) 



Recall that L is called regular when the symmetric matrix [d 2 L/dq % dqi] is every- 
where nonsingular. If L is regular, the Euler-Lagrange equations are second order 
ordinary differential equations for the required curves. 



The standard geometric setting. The action ( |2.lD is independent of the choice 
of coordinates, and thus the Euler-Lagrange equations are coordinate independent 
as well. Consequently, it is natural that the Euler-Lagrange equations may be 
intrinsically expressed using the language of differential geometry. This intrinsic 
development of mechanics is now standard, and can be seen, for example, in Arnold 
[1978], Abraham and Marsden [1978], and Marsden and Ratiu [1994]. 

The canonical 1-form 9q on the 2n-dimensional cotangent bundle of Q, T*Q 
is defined by 

Q (a q )w aq = a q ■ Tir Q w aq , a q £ T*Q, w aq £ T aq T*Q, 

where ttq : T*Q — > Q is the canonical projection. The Lagrangian L intrinsically 
defines a fiber preserving bundle map FL : TQ — > T*Q, the Legendre transfor- 
mation, by vertical differentiation: 



d 

¥L(v q )w q = — 



L(v q + ew q ) 



e=0 
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We define the Lagrange 1-form on TQ, the Lagrangian side, by pull-back 9l = 
¥L*9q, and the Lagrange 2-form by ujl = —dO^. We then seek a vector field Xe 
(called the Lagrange vector field) on TQ such that Xe J u>l = dE, where the 
energy E is defined by E{v q ) = ¥L(v q )v q — L(v q ). 

If FL is a local diffeomorphism then Xe exists and is unique, and its integral 
curves solve the Euler-Lagrange equations (p.4|). In addition, the flow Ft of Xe 
preserves w^; that is, F^uj^ = u^. Such maps are symplectic, and the form ul 
is called a symplectic 2-form. This is an example of a symplectic manifold: a 
pair (M, u) where M is a manifold and uj is closed nondegenerate 2-form. 

Despite the compactness and precision of this differential-geometric approach, 
it is difficult to motivate and, furthermore, is not entirely contained on the La- 
grangian side. The canonical 1-form 9q seems to appear from nowhere, as does 
the Legendre transform FL. Historically, after the Lagrangian picture on TQ was 
constructed, the canonical picture on T*Q emerged through the work of Hamilton, 
but the modern approach described above treats the relation between the Hamil- 
tonian and Lagrangian pictures of mechanics as a mathematical tautology, rather 
than what it is — a discovery of the highest order. 



The variational approach. More and more, one is finding that there are advan- 
tages to staying on the "Lagrangian side". Many examples can be given, but the 
theory of Lagrangian reduction (the Euler-Poincare equations being an instance) is 
one example (see, for example, Marsden and Ratiu [1994] and Holm, Marsden and 
Ratiu [1998a,b]); another, of many, is the direct variational approach to questions 
in black hole dynamics given by Wald [1993]. In such studies, it is the variational 
principle that is the center of attention. 

We next show that one can derive in a natural way the fundamental differential 
geometric structures, including momentum mappings, directly from the variational 
approach. This development begins by removing the boundary condition Sq(a) = 
Sq(b) = from (|2.3| ). Equation ( ^ ) becomes 

J AT \ AT b 

(2.5) 



where the left side now operates on more general dq (this generalization will be 
described in detail in Section (4)), while the last term on the right side does not 
vanish. That last term of ( |2.5| ) is a linear pairing of the function dL/dq 1 , a function 
of q % and q % , with the tangent vector 5q % . Thus, one may consider it to be a 1-form 
on TQ; namely the 1-form (dL/dq l )dq l . This is exactly the Lagrange 1-form, and 
we can turn this into a formal theorem/definition: 



-<k-2 

ping DelL : Q — > T*Q, defined on the second order submanifold 



Theorem 2.1 Given a C k Lagrangian L, k > 2, there exists a unique C k " map- 



2, 



a C 2 curve in Q 
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ofTTQ, and a unique C 1-form 6l on TQ, such that, for all C variations q € (t), 



dS(q(t))-Sq(t) 



where 



Sq(t) = 



de 



D EL L 



£q 
dt 2 



Sq(t) = 



>dt+ e L 



e=0 



de 



e=0 



■Sq 



(2.6) 



t=o 



The 1-form so defined is called the Lagrange I- form. 



and the 



Indeed, uniqueness and local existence follow from the calculation (£ 

coordinate independence of the action, and then global existence is immediate. Here 
then, is the first aspect of our method: 

Using the variational principle, the Lagrange 1-form 0i is the "boundary 
part" of the the functional derivative of the action when the boundary is 
varied. The analogue of the symplectic form is the (negative of) the 
exterior derivative of 6^. 

For the mechanics example being discussed, we imagine a development wherein 9l 
is so defined and we define ujl = —d6i- 

Lagrangian flows are symplectic. One of Lagrange's basic discoveries was that 
the solutions of the Euler-Lagrange equations give rise to a symplectic map. It is a 
curious twist of history that he did this without the machinery of either differential 
forms, of the Hamiltonian formalism or of Hamilton's principle itself. (See Marsden 
and Ratiu [1994] for an account of some of this history.) 

Assuming that L is regular, the variational principle then gives coordinate in- 
dependent second order ordinary differential equations, as we have noted. We tem- 
porarily denote the vector field on TQ so obtained by X, and its flow by Ft. Our 
further development relies on a change of viewpoint: we focus on the restriction of 
S to the subspace Cl of solutions of the variational principle. The space Cl may be 
identified with the initial conditions, elements of TQ, for the flow: to v q G TQ, we 
associate the integral curve s ^> F s (v q ), s £ [0,t]. The value of S on that curve is 
denoted by St, and again called the action. Thus, we define the map St : TQ — > K 
by 



St(Vn 



L(q{s),q(s))ds, 



(2.7) 



where (q{s),q(s)) = F s (v q ). The fundamental equation (|2.6|) becomes 



dS t (v q )w v 



.{Ft{v q )) 



d_ 



F t (v e g ) - 9 L (y q ) -w Vq , 



e=0 



where e i— > v q is an arbitrary curve in TQ such that v q = v q and (d/de)\ov q = w v 
We have thus derived the equation 



ds t = R*e L - e L . 



(2i 
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Taking the exterior derivative of (^]^) yields the fundamental fact that the flow of 
X is symplectic: 

o = dds t = d(F t *e L - e L ) = -f*u l + lj l 

which is equivalent to 

F t *uJ L = uj l . 

This leads to the following: 

Using the variational principle, the fact that the evolution is symplectic 
is a consequence of the equation d 2 = 0, applied to the action restricted 
to the space of solutions of the variational principle. 

In passing, we note that ( |2.8j ) also provides the differential-geometric equations 
for X. Indeed, one time derivative of fl2,8| ) and using ( |2.7D gives dL = £x#c, so that 

X j lo l = -X j d9 L = -£ X L + d(X j Q L ) = d{X j9 L -L) = dE, 

if we define E = X J Ol — L. Thus, we quite naturally find that X = Xe- 

Of course, this set up also leads directly to Hamilton-Jacobi theory, which was 
one of the ways in which symplectic integrators were developed (see McLachlan and 
Scovel [1996] and references therein.) However, we shall not pursue the Hamilton- 
Jacobi aspect of the theory here. 

Momentum maps. Suppose that a Lie group G, with Lie algebra g, acts on Q, 
and hence on curves in Q, in such a way that the action S is invariant. Clearly, G 
leaves the set of solutions of the variational principle invariant, so the action of G 
restricts to Cl, and the group action commutes with Ft. Denoting the infinitesimal 
generator of £ G q on TQ by £tq, we have by Q2.8p , 

o = e TQ j ds t = i TQ j (F t *e L - e L ) = F t * (£ T q j e L ) - e TQ j e L . (2.9) 

For £ 6 g, define Jg : TQ — > M by = £tq -i Ol- Then says that is an 
integral of the flow of Xe- We have arrived at a version of Noether's theorem (rather 
close to the original derivation of Noether): 

Using the variational principle, Noether's theorem results from the in- 
finitesimal invariance of the action restricted to space of solutions of 
the variational principle. The conserved momentum associated to a Lie 
algebra element £ is J§ = £tq -I Ol, where Ol is the Lagrange one-form. 

Reformulation in terms of first variations. We have just seen that symplec- 
ticity of the flow and Noether's theorem result from restricting the action to the 
space of solutions. One tacit assumption is that the space of solutions is a mani- 
fold in some appropriate sense. This is a potential problem, since solution spaces for 
field theories are known to have singularities (see, e.g., Arms, Marsden and Moncrief 



8 



[1982]). More seriously there is the problem of finding a multisymplectic analogue 
of the statement that the Lagrangian flow map is symplectic, since for multisym- 
plectic field theory one obtains an evolution picture only after splitting spacetime 
into space and time and adopting the "function space" point of view. Having the 
general formalism depend either on a spacetime split or an analysis of the associated 
Cauchy problem would be contrary to the general thrust of this article. We now 
give a formal argument, in the context of Lagrangian mechanics, which shows how 
both these problems can be simultaneously avoided. 

Given a solution q(t) £ Cl, a first variation at q(t) is a vector field V on Q 
such that t i — > FY ° q(t) is also a solution curve (i.e. a curve in Cl)- We think of 
the solution space Cl as being a (possibly) singular subset of the smooth space of 
all putative curves C in TQ, and the restriction of V to q(t) as being the derivative 
of some curve in Cl at q(t). When Cl is a manifold, a first variation is a vector at 
q(t) tangent to Cl- Temporarily define a = dS — 6l where by abuse of notation Ol 
is the one form on C defined by 

9 L {q{t))5q{t) = 9 L (b)5q(b) - 9 L (a)6q(a). 

Then Cl is defined by a = and we have the equation 

dS = a + 9 L , 

so if V and W are first variations at q(t), we obtain 

= V AW Ad 2 S = V AW Ada + V AW Ad9 L - (2.10) 

We have the identity 

da(V,W)(q(t)) = V(a(W)) - W(a(V)) - a([V,W]), (2.11) 

which we will use to evaluate ( |2.10| ) at the curve q(t). Let FY denote the flow of V, 
define qY(t) = FY (q(t)), and make similar definitions with W replacing V. For the 
first term of ( |2.11| ), we have 



V(a{W))[q(t)) = ± 



a(W)(qY), 

e=0 



which vanishes, since a is zero along qY for every e. Similarly the second term 



of (2.11) at q(t) also vanishes, while the third term vanishes since a (<?(£)) = 0. 



Consequently, symplecticity of the Lagrangian flow Ft may be written 

V A WAd6 L = 0, 

for all first variations V and W . This formulation is valid whether or not the solution 
space is a manifold, and it does not explicitly refer to any temporal notion. Similarly, 
Noether's theorem may be written in this way. Summarizing, 
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Using the variational principle, the analogue of the evolution is symplec- 
tic is the equation d 2 S = restricted to first variations of the space of 
solutions of the variational principle. The analogue of Noether's theorem 
is infinitesimal invariance of dS restricted to first variations of the space 
of solutions of the variational principle. 

The variational route to the differential-geometric formalism has obvious ped- 
agogical advantages. More than that, however, it systematizes searching for the 
corresponding formalism in other contexts. We shall in the next sections show how 
this works in the context of discrete mechanics, classical field theory and multisym- 
plectic geometry. 



3 Veselov Discretizations of Mechanics 

The discrete Lagrangian formalism in Veselov [1988], [1991] fits nicely into our vari- 
ational framework. Veselov uses QxQ for the discrete version of the tangent bundle 
of a configuration space Q; heuristically, given some a priori choice of time interval 
At, a point (gi, go) £ Q x Q corresponds to the tangent vector (qi — q )/At. Define 
a discrete Lagrangian to be a smooth map L : Q x Q = {gi,go} — ► M, and the 
corresponding action to be 

n 

S = ^L(g fe ,9fc-i)- (3.1) 
k=i 

The variational principle is to extremize S for variations holding the endpoints go 
and q n fixed. This variational principle determines a "discrete flow" F : Q x Q — > 
QxQ by F(qi,qo) = (g2, gi), where g2 is found from the discrete Euler- Lagrange 
equations (DEL equations): 

dL . . dL , . 

t- gi,9o ) + -3- 92,91 =0. 3.2 
dqi dq 

In this section we work out the basic differential-geometric objects of this discrete 
mechanics directly from the variational point of view, consistent with our philosophy 
in the last section. 

A mathematically significant aspect of this theory is how it relates to integrable 
systems, a point taken up by Moser and Veselov [1991]. We will not explore this 
aspect in any detail in this paper, although later, we will briefly discuss the reduc- 
tion process and we shall test an integrator for an integrable pde, the sine-Gordon 
equation. 

The Lagrange 1-form. We begin by calculating dS for variations that do not fix 
the endpoints: 

dS(q ,--- ,9„) • (<%>,••• ,Sq n ) 
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^2 (J^(qk+i,qk)Sqk+i + ^(qk+i,qk)&qi 
Q—(qk,qk-i)Sqk + Q—(qk+i,qk)5qk 

fe=l 91 fc=0 ?0 



k=l 

dL dL 

+ -^—(qi,qo)$qo + -K—(qn,qn-i)8qn- (3.3) 

dq dqi 

It is the last two terms that arise from the boundary variations (i.e. these are the 
ones that are zero if the boundary is fixed), and so these are the terms amongst 
which we expect to find the discrete analogue of the Lagrange 1-form. Actually, 
interpretation of the boundary terms gives the two 1-forms on Q x Q 

dL 

9 L (qi,qo) • (Sqi,Sq ) = j—(qi,qo)Sqo, (3.4) 

oq 

and 

dL 

6l(qi,qo) ■ (Sqi,Sq ) = j—(qi,qo)8qi, (3.5) 

oqi 

and we regard the pair {6~ ,0 + ) as being the analogue of the one form in this situa- 
tion. 

Symplecticity of the flow. We parameterize the solutions of the variational 
principle by the initial conditions (qi,qo), and restrict S to that solution space. 
Then Equation ( |3.3| ) becomes 

dS = 91 +F*6+. (3.6) 

We should be able to obtain the symplecticity of F by determining what the equation 
ddS = means for the right-hand-side of Q3.6j). At first, this does not appear to 
work, since ddS = gives 

F*(dB+) = -Mi, (3.7) 

which apparently says that F pulls a certain 2-form back to a different 2-form. The 
situation is aided by the observation that, from ( |3.4| ) and ( |3.5| ), 

9i + 9+ = dL, (3.8) 

and consequently, 

d6 L + del = 0. (3.9) 
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Thus, there are two generally distinct 1-forms, but (up to sign) only one 2-form. If 
we make the definition 



UL = dei = -del, 

then becomes F*ujl = uj^. Equation (3^), in coordinates, gives 



Wi = a in J d % Ad( fv 

which agrees with the discrete symplectic form found in Veselov [1988], [1991]. 

Noether's Theorem. Suppose a Lie group G with Lie algebra q acts on Q, and 
hence diagonally on Q x Q, and that L is G-invariant. Clearly, S is also G-invariant 
and G sends critical points of S to themselves. Thus, the action of G restricts to 
the space of solutions, the map F is G-equivariant, and from (|3.6|), 



= £q x q JdS = £q x q J 0- L + £qxq J (F*9+), 
for £ 6 q, or equivalently, using the equivariance of F, 

H QxQ jei = -F*(Z QxQ J0 + ). (3.10) 



Since L is G-invariant, ( plf gives £qxQ J 9 L = —£,q x q -I which in turn con- 
verts ( |3.10| ) to the conservation equation 

£ Qx qJ0+ = F*(£ QxQ j0+). (3.11) 

Defining the discrete momentum to be 

we see that ( |3 . 1 1| ) becomes conservation of momentum. A virtually identical deriva- 
tion of this discrete Noether theorem is found in Marsden and Wendlant [1997] . 

Reduction. As we mentioned above, this formalism lends itself to a discrete ver- 
sion of the theory of Lagrangian reduction (see Marsden and Scheurle [1993a, b], 
Holm, Marsden and Ratiu [1998a] and Cendra, Marsden and Ratiu [1998]). This 
theory is not the focus of this article, so we shall defer a brief discussion of it until 
the conclusions. 



4 Variational Principles for Classical Field Theory 

Multisymplectic geometry. We now review some aspects of multisymplectic 
geometry, following Gotay, Isenberg and Marsden [1997] and Marsden and Shkoller 
[1997]. 

We let ttxy ■ y —> X be a fiber bundle over an oriented manifold X. Denote 
the first jet bundle over Y by J 1 (Y) or J X Y and identify it with the affine bundle 
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over Y whose fiber over y G Y x := 7r XY (x) consists of AS(T x X,T y Y), those linear 
mappings 7 : T X X — ► T y Y satisfying 

T7Txy 07 = Identity on T X X. 

We let dimA = n + 1 and the fiber dimension of Y be TV. Coordinates on X are 
denoted /U = 1, 2, . . . , n, 0, and fiber coordinates on Y are denoted by y" 4 , A = 
1, . . . , N. These induce coordinates v A ^ on the fibers of J 1 (y). If 4> : X — > Y is a 
section of 7Txr, its tangent map at 1 £ I, denoted T x cj), is an element of J l {Y)^ x y 
Thus, the map x ^ T x (f> defines a section of J 1 (Y) regarded as a bundle over X. 
This section is denoted j 1 (4>) or j 1 ^ and is called the first jet of <f>. In coordinates, 
j 1 {4>) is given by 

x ^ (a^^&^V)), (4.1) 

where = d/dx u . 

Higher order jet bundles of Y, J m (Y), then follow as J x (- • •(J 1 (Y)). Analogous 
to the tangent map of the projection t^y,j 1 (y)^ '^' k y,j 1 (y) '■ TJ l (Y) — ► TY, we may 
define the jet map of this projection which takes J 2 (Y) onto J 1 (Y) 

Definition 4.1 Xei 7 G «/ 1 (Y) so that ^x,J 1 (Y)(l) = x - Then 

j7r YJ1(Y) : AfT^A, T^Y)) - Aff(T x X, T7r YiJl(Y) • T^Y)). 

VFe define the subbundle Y" of J 2 (Y) over X which consists of second- order jets so 
that on each fiber 

Y; = { S G J 2 (Y) 7 I J7r yiJl(y) ( S )= 7 }. 

In coordinates, if 7 G «/ 1 (Y) is given by (x^,y A ,v A ^), and s G J 2 (Y) 7 is given 
by (x^,y A ,v A fl ,w A ^,K A fll/ ), then s is a second-order jet if v A M = Thus, 
the second jet of 4> G r(7rxy), j 2 {<fi), given in coordinates by the map x' 1 i-> 
(x^, t^" 4 , d,y(f) A , d/j,d u (f) A ), is an example of a second-order jet. 

Definition 4.2 The dual jet bundle J X (Y)* is the vector bundle over Y whose 
fiber at y G Y x is the set of affine maps from J 1 (Y) y to A n+1 (X) x , the bundle of 
(n + l)-forms on X. A smooth section of J 1 (Y)* is therefore an affine bundle map 
o/J^Y) to A n+1 (X) covering ir xy ■ 

Fiber coordinates on J 1 (Y)* are (p,PA^), which correspond to the affine map given 
in coordinates by 

A- {p + PA»v A ^d n+1 x, (4.2) 

where d n+1 x = dx 1 A • • • A dx 11 A dx°. 

Analogous to the canonical one- and two-forms on a cotangent bundle, there exist 
canonical (n + 1)- and (n + 2)-forms on the dual jet bundle J X (Y)*. In coordinates, 
with d n x^ := <9 M J d n+1 x, these forms are given by 

6 = PA tl dy A A d n Xfl + pd n+1 x and Q = dy A A dp A ^ A d n x^ - dp A d n+1 x. (4.3) 
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A Lagrangian density £ : ^{Y) — ► A n+1 (A) is a smooth bundle map over X. 
In coordinates, we write 



£( 1 ) = L(x»,y A ,v A fl )d n+1 x. 



(4.4) 



The corresponding covariant Legendre transformation for £ is a fiber preserving 
map over Y, F£ : J l (Y) — ► J 1 (y)*, expressed intrinsically as the first order vertical 
Taylor approximation to £: 



F£( 7 ) • i = £(7) + 



de 



£(7 + 5(7' -7)) 



(4.5) 



e=0 



where 7,7' € J 1 (y) J/ . A straightforward calculation shows that the covariant Leg- 
endre transformation is given in coordinates by 



PA h 



dL 

dv A „ 



and p = L — 



(4.6) 



We can then define the Cartan form as the (n + l)-form 0£ on J 1 (l") given 



by 

e £ = (F£)*e, 

and the (n + 2)-form Qc by 

f) £ = -d® c = (¥C)*n, 
with local coordinate expressions 



(4.7) 
(4.8) 



dL 

dv A , 



dy A A d n Xn + [L - 



dL 

dv A ,, 



v A » ) d n+1 x, 



Q c = dy A A d 



dL 

dv A , 



A d n x fl - d 



L 



dL A ' 

dv\ v r 



(4.9) 



A d n+1 



x. 



This is the differential-geometric formulation of the multisymplectic structure. 
Subsequently, we shall show how we may obtain this structure directly from the 
variational principle, staying entirely on the Lagrangian side J x (y). 



The multisymplectic form formula. In this subsection we prove a formula that 
is the multisymplectic counterpart to the fact that in finite-dimensional mechanics, 
the flow of a mechanical system consists of symplectic maps. Again, we do this by 
studying the action function. 

Definition 4.3 Let U be a smooth manifold with (piecewise) smooth closed bound- 
ary. We define the set of smooth maps 

C°° ={(/>: U — > Y I ttxy : U — > X is an embedding}. 

For each 4> G C°° , we set <px '■= kxy°4> an d Ux '■= ^xy °4>(U) so that (fix '■ U — > Ux 
is a diffeomorphism. 
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We may then define the infinite-dimensional manifold C to be the closure of C°° 
in either a Hilbert space or Banach space topology. For example, the manifold C 
may be given the topology of a Hilbert manifold of bundle mappings, H S (U, Y), (U 
considered a bundle with fiber a point) for any integer s > (n + l)/2, so that the 
Hilbert sections (f> o <f) x x in Y are those whose distributional derivatives up to order s 
are square-integrable in any chart. With our condition on s, the Sobolev embedding 
theorem makes such mappings well defined. Alternately, one may wish to consider 
the Banach manifold C as the closure of C°° in the usual C fc -norm, or more generally, 
in a Holder space C fc+Q -norm. See Palais [1968] and Ebin and Marsden [1970] for 
a detailed account of manifolds of mappings. The choice of topology for C will not 
play a crucial role in this paper. 

Definition 4.4 Let Q be the Lie group of ir xy -bundle automorphisms t]y covering 
diffeomorphisms nx, with Lie algebra q. We define the action <3? : Q x C — > C by 

Hvy, <f>) = vy © ((f) o cj) x l ) o n' 1 . 

Furthermore, if <f) o (f)^ 1 G T(ttu x ,y), then $(r/y,0) G ^(^rj X (u x ),y)' The tangent 
space to the manifold C at a point 4> is the set T^C defined by 

{V G C°°(X, TY) | it y ,ty °V = 4>, &Tn X Y o V = V x , a vector field on X} . 

(4.10) 

Of course, when these objects are topologized as we have described, the definition 
of the tangent space becomes a theorem, but as we have mentioned, this functional 
analytic aspect plays a minor role in what follows. 

Given vectors V, W G T^C we may extend them to vector fields V, W on C by 
fixing vector fields v, w G TY such that V = v o ((f> o (f) x l ) and W = w o ((f) o (f)^ 1 ), 
and letting V p = v o (p o p x l ) and W p = w o (po p^ 1 ). Thus, the flow of V on C 
is given by ^(r] Y ,p) where n Y covering rj x is the flow of v. The definition of the 
bracket of vector fields using their flows, then shows that 

[V,W](p) = [v,w]o(po Px 1 ). 
Whenever it is contextually clear, we shall, for convenience, write V for v o(<p o (f)^ 1 ). 
Definition 4.5 The action function S on C is defined as follows: 

S(4>)= l £(j 1 (0o0 x 1 ))/ O ra/^GC. (4.11) 
JU X 

Let A i— > r) Y be an arbitrary smooth path in Q such that rfy = e, and let V G T^C 
be given by 



V = — 

dX 



A=0 « A 



Vx o 4>- (4-12) 

A=0 
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Definition 4.6 We say that <fr is a stationary point, critical point, or ex- 
tremum of S if 



d_ 
dX 



S($(r£,0)) = O. 



A=0 



Then, 



dSs-V 



_d_ 

dX 
d_ 
dX 



A=0 J Vx 0< t>x(U) 



J 



(4.13) 



(4.14) 



- 1 )*i 1 (^)*9 £ , 



\=0J<f>x(U) 

where we have used the fact that C(z) = z*Qc f° r an holonomic sections z of J l (Y) 
(see Corollary \fl below), and that 

j x (n Y o (j) o (j)' 1 o n^}) = j x {tiy) o j x {(j) o 0" 1 ) o jy" 1 . 
Using the Cartan formula, we obtain that 



-^*\j\v)jn c ) 



J x 



+ 



J 



dUx 



y[j L (v)je c ]. 



(4.15) 



Hence, a necessary condition for G C to be an extremum of S is that the first 
term in ( |4.15j ) vanish. One may readily verify that the integrand of the first term in 
( 4,15| ) is equal to zero whenever j (V) is replaced by W G T^iY) which is either 
tt^jw) -vertical or tangent to j 1 ^ o (p^) (see Marsden and Shkoller [1997]), so that 
using a standard argument from the calculus of variations, j x (<p o (^• 1 )*[Wjfi£] must 
vanish for all vectors W on ^(Y) in order for <p to be an extremum of the action. 
We shall call such elements (j) £ C covering (j>x, solutions of the Euler-Lagrange 
equations. 



Definition 4.7 We let 

V = {(/) G C | j'^o^r^Jflr] = for all W G TJ l (Y)}. 
In coordinates, ((f) o 4>^) A is an element ofV if 



(4.16) 



dL 

dy A 



<9a^ \dv A 



We are now ready to prove the multisymplectic form formula, a generalization 
of the symplectic flow theorem, but we first make the following remark. If V is 
a submanifold of C, then for any <p G V, we may identify T^P with the set {V G 
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T^C | j l ((f> o <!>](')* £ji(y) [W -I Sic] = for all W G TJ 1 ^)} since such vectors 
arise by differentiating ^ |e=oJ 1 (0 e ° ^x -1 )*!^- 1 ^c] = 0, where <^ e is a smooth 
curve of solutions of the Euler-Lagrange equations in V (when such solutions exist). 
More generally, we do not require V to be a submanifold in order to define the first 
variation solution of the Euler-Lagrange equations. 

Definition 4.8 For any eft € V ,we define the set 

T = {V£T^C \ j\<Po0 x 1 )*£ jl{v) [WjSl c } = O for all W E TJ X (Y)}. (4.17) 

Elements of T solve the first variation equations of the Euler-Lagrange equations. 



Theorem 4.1 (Multisymplectic form formula) If (ft £ V, then for all V and 
W in T, 



j 1 ((fto ( ft x 1 y[j 1 (v)jj 1 (w)jnc} = o- 



(4.18) 



Proof. We define the 1-forms ol\ and a 2 on C by 

ax((ft)-V:= [ ^((fto^yif^jSlc] 
Ju x 



and 



a 2 ((ft)-V :-- 



dU x 



j\(ftocft x U *' V 



\jHv)j®c], 



so that by ( |4.15[ ), 

dStj, ■ V = ai{(ft) ■ V + a 2 (<A) • V for all V € T^C. 
Recall that for any 1-form a on C and vector fields V, W on C, 
da(V, W) = V[a{W)] - W[a{V)] - a([V, W]). 



(4.19) 



(4.20) 



We let (ft e = r/y o (ft be a curve in C through (ft, where n Y is a curve in Q through the 
identity such that 

W = ^-l=or]Y and W e J 7 , 
de 



and consider equation ( 4.19 ) restricted to all V € T . 
Thus, 



d{a 2 {V)){cp)-W = j e 

d_ 
~de 



MV){(ft e )) 



e=0 



=0 JdUy, 



jHcfto^yfiny^^jQc] 
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JdUx 

- [ j 1 (^ocf ) ^r\j 1 (W)jd(j\V)jQc)} 

JdUx 

+ f j 1 (4>o^ 1 rd\j 1 (w)jj 1 (v)je c }, 

JdUx 

where the last equality was obtained using Cartan's formula. Using Stoke's theorem, 
noting that ddU is empty, and applying Cartan's formula once again, we obtain that 

d(a 2 (cf>)(V))-W = [ j 1 {ct>oct>- 1 y\j\W)Aj\V)An c \ 
JdUx 

- [ /(^^rb'W^wed, 

JdU x 

and 

d(a 2 (4>)(W))-V = [ j\<l ) o ( p x 1 y[j\v)jj\W)jQc} 

JdU x 

JdUx 

Also, since [j 1 {V) , j 1 {W)\ = jH^W}), we have 

c*2((f>W,w}) = [ j 1 Wo^ 1 rb' 1 (v)>j' 1 W]je £ . 

JdUx 

Now 

\j\v),j\w)\ a@ c = z^tiHw) j @c) - j\w) j %i (y) e £ , 

so that 

dn 2 {<t>){V,W) = 2 I j 1 ( ( f>o ( t>- 1 )*[j 1 (V)jj\W)jnc} 
JdUx 

+ [ j\<t> o <t>xT\j\V) j z jHw) e c - Z^tiHW) J &c)}. 

JdUx 

But 

£j\v)(i\W) J @c) = d{j l {V)Aj\W) J B c ) +j\V) J d(j\W) J @ c ) 

and 

j\v) j %i (w) e £ = j\v) j dO' 1 ^) J ©/:) +j\v) J i 1 W J 

Hence, 

JdU x 
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/ j i (<t>o^ i rd(j i (v)jj i (w)j& c ). 

JdU x 



The last term once again vanishes by Stokes theorem together with the fact that 
ddU is empty, and we obtain that 

da 2 (<t>)(v,w) = 3 / j 1 (4>o4 >x 1 y(j 1 (v)jj 1 (w)jn c ). (4.21) 

JdUx 

We now use (4.20) on a\. A similar computation as above yields 

d(a 1 ( ( f>)-v)-w= [ jHfo&yzpwtfQO-inc] 

JUx 

which vanishes for all (f) G V and W G T . Similarly, d(ai(4>) ■ W) ■ v = for all 
4> G V and V G ^. Finally, ai(» = for all <f> G P. 
Hence, since 

= ddS(<£) ( V, WO = dai (0) ( V, WO + rf «2 (</>)(V, W), 
we obtain the formula ( |4.18| ), ■ 



Symplecticity revisited. Let S be a compact oriented connected boundaryless 
n-manifold which we think of as our reference Cauchy surface, and consider the space 
of embeddings of S into X, Emb(S,X); again, although it is unnecessary for this 
paper, we may topologize Emb(S,X) by completing the space in the appropriate 
C k or .fP-norm. 

Let B be an m-dimensional manifold. For any fiber bundle ttbk '■ K — > B, we 
shall, in addition to T(itbk), use the corresponding script letter K, to denote the 
space of sections of ttbk- The space of sections of a fiber bundle is an infinite- 
dimensional manifold; in fact, it can be precisely defined and topologized as the 
manifold C of the previous section, where the diffeomorphisms on the base manifold 
are taken to be the identity map, so that the tangent space to K at a is given simply 
by 

T a K, = {W : B -^VK \k k ,TK oW = a}, 

where VK denotes the vertical tangent bundle of K. We let ^K,L(VK,k m (B)) '■ 
L(VK,A m (B)) — > K be the vector bundle over K whose fiber at k G K x , x = 
^BK(k), is the set of linear mappings from V^K to A m (B) x . Then the cotangent 
space to K, at a is defined as 

T*/C = {tt : B - L(VK, A m (B)) \ n KMVK:An+ i (B)) 07T = a}. 

Integration provides the natural pairing of T*/C with T a K,\ 

(TT,V)= [ TT- V. 

Jb 

In practice, the manifold B will either be X or some (n + l)-dimensional subset of 
X, or the n-dimensional manifold S T , where for each r G Emb(E,JT), S T := r(S). 
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We shall use the notation Y T for the bundle ttt, t ,y, and y T for sections of this bundle. 
For the remainder of this section, we shall set the manifold C introduced earlier to 

y. 

The infinite-dimensional manifold y r is called the r- configuration space, its 
tangent bundle is called the r-tangent space, and its cotangent bundle T*y r is 
called the r -phase space. Just as we described in Section |2|, the cotangent bundle 
has a canonical 1-form 8 T and a canonical 2-form uj t . These differential forms are 
given by 

6 T (<p,n) • V = [ iT(Tiry T . T *y T ■ V) and uj t = -d6 T , (4.22) 

where (<p,Tr) G y r , V G T^^T*y r , and ^y T ,T*y T ■ T*y T — > >V is the cotangent 
bundle projection map. 

An infinitesimal slicing of the bundle ttxy consists of Y T together with a vector 
field £ which is everywhere transverse to Y T , and covers £x which is everywhere 
transverse to S T . The existence of an infinitesimal slicing allows us to invariantly 
decompose the temporal from the spatial derivatives of the fields. Let <fi G y, 
ip := </>|s r , and let i T : T, T — » X be the inclusion map. Then we may define the map 
P c taking j 1 (y) T to j l {y> T ) X r(7r Er) vy r ) over ^ T by 

PctiHt) ° V) = (i 1 ^)^) where <p := £ c 0. (4.23) 

In our notation, j 1 (y) T is the collection of restrictions of holonomic sections of 
J l (Y) to E r , while j (y T ) are the holonomic sections of 7r STj ji(y). It is easy to see 
that /?£ is an isomorphism; it then follows that (5^ is an isomorphism of j 1 (y) T with 
Ty T , since j (tp) is completely determined by cp. This bundle map is called the jet 
decomposition map, and its inverse is called the jet reconstruction map. Using this 
map, we can define the instantaneous Lagrangian. 

Definition 4.9 The instantaneous Lagrangian L T ^ : Ty r — ► R is given by 

L T)( {<p, ( p)= f e T [Q x AC(^\j 1 ^),^)} (4.24) 

for all (ip, ip) G Ty r . 

The instantaneous Lagrangian L T £ has an instantaneous Legendre transform 

¥L T>C : Ty T ^ T*y T ; (ip, <p) ^ vr) 

which is defined in the usual way by vertical fiber differentiation of L T ^ (see, for 
example, Abraham and Marsden [1978]). Using the instantaneous Legendre trans- 
formation, we can pull-back the canonical 1- and 2-forms on T*y r . 

Definition 4.10 Denote, respectively, the instantaneous Lagrange 1- and 2-forms 
on Ty r by 

0^ = ¥L* x e T and ^ = -d9%. (4.25) 
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Alternatively, we may define 6^ using Theorem 2.1, in which case no reference to 
the cotangent bundle is necessary. 

We will show that our covariant multisymplectic form formula can be used 
to recover the fact that the flow of the Euler-Lagrange equations in the bundle 
7r Emb(s,x),u TeEmb(E X )Ty T is symplectic with respect to u>^. To do so, we must relate 
the multisymplectic Cartan (n + 2)-form f^£ on J 1 (1 A ) with the symplectic 2-form 
ujL on Ty T . 

Theorem 4.2 Let be the canonical 1-form on j 1 (y) T given by 

ef(j 1 (^)o V ).y= / i^^rivjQc), (4.26) 

where j\<j>) o i T G j l {y) T , V G T^^J 1 ^. 

(a) If the 2-form on j 1 (y) T is defined by fif := — dQ^ , then for V, W G 



/ i* T j\<j>)*[W AV ASlc\ (4.27) 



(b)Let the diffeomorphism sx '■ S x R — > X be a slicing of X such that for A 6 t, 

Sa := sx(S x {A}) and S A := r A (S), 

where t\ G Emb(S,X) is given 6y ta(x) = sx(x,X). For any <p £ V, let V, W G 
T^nf so that for each r G Emb(S,A") ; j 1 V T ,j 1 W T G 2)i ( ( p) i r j 1 (y)r, and let 
T Ai) T A 2 ^ Emb(E, X). Then 

^(i 1 ^, j 1 ^) = 0^(j% 2 , /W^). (4.28) 

Proof. Part (a) follows from the Cartan formula together with Stokes theorem 
using an argument like that in the proof of Theorem |4.1| , 

For part (b), we recall that the multisymplectic form formula on y states that 
for any subset Ux C X with smooth closed boundary and vectors V, W G T^y n J 7 , 

[ jH4>r[jHv)jj\w)jn c ] = o. (4.29) 

JdUx 

Let 

Ux = UA6[Ax,Aa] S A- 
Then df/x = '^x 1 — ^a 2 , so that ( 4.29| ) can be written as 

- / /(0oi j/^ TAi J0 £ ] 

= ^ (i% a , j'w TXi ) - (i x K A2 ,i 1 H^ TA2 ), 

which proves ( 4.28| ). ■ 
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Theorem 4.3 The identity 6^ = fiffi holds. 

Proof. Let W £ Tjir^^j (y) T , which we identify with w o <p o i T1 where w is a 
ttx,ji(y) -vertical vector. Choose a coordinate chart which is adapted to the slicing 
so that d \Y T = C With w = (0, W A , W A ), we see that 



Now, from (4.24) we get 



%M) = d ^y A 

d 



i*[d AL(x»A A A A Jd n+l x®dy A ] 

dL 



B ±B \i,A 



T 3u A 



,(/>" tlJ ,)dy A ®d n x , 



where we arrived at the last equality using the fact that y A = v A q in this adapted 
chart. Since (T/? c • W) A = W A , we see that Q$ ■ W = 6^ ■ (T(3 ( ■ W), and this 
completes the proof. ■ 

Let the instantaneous energy E T ^ associated with L T ^ be given by 

E T)( (<p, <p) = FL T)f (0) • tp - L T>c (<p, <j>), (4.30) 
and define the "time" -dependent Lagrangian vector field Xe t( by 

X Eti( =dE T>c . 

Since U Tg Emb(s,x)^3 ; T over Emb(£, X) is infinite-dimensional and is only weakly 
nondegenerate, the second-order vector field Xe t( does not, in general, exist. In 
the case that it does, we obtain the following result. 

Corollary 4.1 Assume Xe t( exists and let F T be its semiflow, defined on some 
subsetV of the bundle U rG Emb(s,Js:)^3 ; r over Emb(T,, X). Fixf so that Ff(ipi,<pi) = 
(^2,^2) where (<pi,ipi) G Ty ri and (^2,^2) G Ty T2 . Then F?u^ = 



Proof. This follows immediately from Theorem |4.2| (b) and Theorem 4.3 and the 
fact that /?£ induces an isomorphism between j 1 (y) T and Ty r . ■ 

Example: nonlinear wave equation. To illustrate the geometry that we have 
developed, let us consider the scalar nonlinear wave equation given by 

^ _ A - AT'(0) = 0, 0Gr(7r X y), (4.31) 



dx°< 
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where A is the Laplace-Beltrami operator and N is a real-valued C°° function of 
one variable. For concreteness, fix n=l so that the spacetime manifold X := IR 2 , 
the configuration bundle Y := 7r K 2 k , and the first jet bundle := 7r K 2 iR 3. 

Equation ( 4.3l| ) is governed by the Lagrangian density 



d^_ 2 
dx° 



dx 1 



+ N(<f>) } dx 1 A dx Q 



(4.32) 



Using coordinates (x , x 1 , c/>, <fi t Q, (p^) for J l {Y), we write the multisymplectic 3-form 
for this nonlinear wave equation on R 2 in coordinates as 



> A d(j) fi A dx 1 — d<f> A #,i A dx° - N'((j))d(j) A dx 1 A dx° 
+4> fi d(f)fi A dx 1 A dx° - A dx 1 A a!x ; (4.33) 



a short computation verifies that solutions of ( [4.31 ) are elements of V, or that 
j\<j> o <f>£)*\W J n c ] = for all W £ TJ l {Y) (see Marsden and Shkoller [1997]). 

We will use this example to demonstrate that our multisymplectic form formula 
generalizes the notion of symplecticity given by Bridges [1997]. Since the Lagrangian 
( 4.32j ) does not explicitly depend on time, it is convenient to identify sections of Y 
as mappings from IR 2 into IR, and similarly, sections of J 1 (Y) as mappings from IR 2 
into IR 3 . Thus, for <j) G r(vrxy), j 1 (4>)(x fl ) := (0(a^), ^o(^), G and if 

we set p M : 



4>> 



then (^3l|) can be reformulated to 



f oiV,o + Jii 




,1 : = 
















10" 




" 4> ' 




" 





-1 " 




" 4> ' 




r N'w 1 


-100 




p° 


+ 













p° 




-p° 


000 




pi 


,0 


1 










pi 


,1 


pi 



(4.34) 



To each degenerate matrix J^, we associate the contact form on IR 3 given by 
0^(^1,1^2) = (J M ni,U2), where U\,U2 £ IR 3 and (•,•) is the standard inner product 
on M 3 . Bridges obtains the following conservation of symplecticity: 

^0 hVW-o),^,!))] + kC^fto),^^)] = 0. (4-35) 
This result is interesting, but has somewhat limited scope in that the vector 



fields in ( 4.35 ) upon which the contact forms act are not general solutions to the 
first variation equations; rather, they are the specific first variation solutions </> iAt . 
Bridges obtains this result by crucially relying on the multi-Hamiltonian structure 
of ( [4.31| ); in particular, the vector (N'(<p), — p Q ,p l ) on the right-hand-side of ( |4.34 ) 
is the gradient of a smooth multi-Hamiltonian function H(<p,p° ,p l ) (although the 
multi-Hamiltonian formalism is not important for this article, we refer the reader to 
Marsden and Shkoller [1997] for the Hamiltonian version of our covariant framework, 
and to Bridges [1997]). Using equation (|4.34 ), it is clear that 



H, 



«> UH<l>,o),jH<l>,i)) and 

so that ( 4.35] ) follows from the relation -ff 0,1 = 



= -u,\j 1 (^ )J 1 (<P,i)) 



23 



Proposition 4.1 The multisymplectic form formula is an intrinsic generalization 



of the conservation law fy.Sdj ); namely, for any V,W G T that are ir x . 



JHY) 



■vertical, 



A [u:\j\Y) ,j\W))] + A [JtfQOjHW))] = 0. (4.36) 

Proof. Let j (V) and j (W) have the coordinate expressions (V,V° ,V X ) and 
(W, W°, W 1 ), respectively. Using ( [4.33 ), we compute 



j l {W) J j x {V) Jn c = {VW° - V°W) dx + (VW 1 - V X W) alt, 
so that with Theorem |4.1| and the definition of w^, we have, for Ux C X, 

u J °(j 1 (V),j\W))dx-u; 1 (j 1 (V),j 1 (W))dt = 0, 



dU x 

and hence by Green's theorem, 

{^o W{v), 3 Hw))} + A [ w W), i 1 (WO)] } ^ a dx° = o. 

Since Ux is arbitrary, we obtain the desired result. ■ 

In general, when V is 7rxy-vertical, .7 (V) has the coordinate expression (V, + 
dV/d(j) ■ 4>,n , but for the special case that V = <f) ;fl , j yt>n) = C/V),/^ and Proposition 
Tl gives 

d d 
O^o [0o0,o,i - <h<f>,o,o\ ~ g^l [0o0,i,i - 010,O,l] = 0, 

which simplifies to the trivial statement that 

0, o iV(0),i - 0,iiV(0),o = 0. 

The variational route to the Cartan form. We may alternatively define the 
Cartan form by beginning with equation ( [4.14| ). Using the infinitesimal generators 
defined in (4.12), we obtain that 



dSck ■ V 



d 




dX 


A=0 


Jv x (Ux) 




d 


1 

Ju x 


dX 



d 



Ctf{*{rk,4>))) 



A=0 



A=0 

+ / £ Vx [£(j 1 (0o0 x 1 ))] 

'Ux 



(4.37) 



Using the natural splitting of TY, any vector V £ T^C may decomposed as 

V = V h + V v , where V h = T(0 o X X ) • V x and V v = V - V h , (4.38) 
where we recall that Vx = Tttxy • V. 
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Lemma 4.1 For any V E T^C, 
dS^ ■ V h -- 

and 

dS* ■ V v -- 



U X 



d\ 



(4.39) 



(4.40) 



A=0 



Proof. The equality ( [4.40 ) is obvious, since the second term in flOg) clearly 
vanishes for all vertical vectors. For vectors V , the first term in (|4.37|) vanishes; 
indeed, using the chain rule, we need only compute that 



d 

dX 



V h -T{cf>o4>- 1 ).Vx 



A=0 



which is zero by ( 4.38 ). We then apply the Cartan formula to the second term in 
( 4.37] ) and note that dC is an (n + 2)-form on the (n + l)-dimensional manifold Ux 
so that we obtain ( 4.39j ). ■ 

Theorem 4.4 Given a smooth Lagrangian density £ : ^(Y) — > A n+1 (X), there 
exist a unique smooth section Del£ S C°°(Y", A n+1 (X) ® T*Y)) and a unique 
differential form 0£ G A n+1 (J 1 (Y)) such that for any V G T^C, and any open 
subset Ux such that [/jfl dX = 0, 

dS 4> -V=( D EL £(j 2 (<po^ x l )).V+ [ j^o^WOOJ®/:]- Ml) 

JUx JdUx 



Furthermore, 

D EL C(j 2 (^ o ■ V = j\<j> o ^ytfiV) J Sl c ] in U x . (4.42) 
In coordinates, the action of the Euler- Lagrange derivative Del£ on Y" is given by 



d 2 L 



dL 



QyBQyA^ 

dv B ,,dv A , 



aUH'P^x 1 )) 



8 2 L 



dx^dv A , 



-(/OW* 1 )) 



-(i 1 (0o^ 1 )).(^o0 x 1 ^ 



1\B 

flU 



dy A A d n+l x 



(4.43) 



while the form ©£ matches the definition of the Cartan form given in and has 
the coordinate expression 



<>L dy A A d n x (i + ( L 



dv A 



dL r A |J" + J T 
dv\ V M ' 



(4.44) 
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Proof. Choose Ux '■= <ftx(U) small enough so that it is contained in a coordinate 
chart, say 0\. In these coordinates, let V = (V\ V 71 ) so that along (^o^" 1 , our 
decomposition fl4.3Sj ) may be written as 



v A - V 1 



d(4> o (f> 



-1\A 



X 



Q x fj, 



d 

dy / 



and equation (|4.40|) gives 



dL . x 



v\Al 



dx^ 



d n+1 x, 
(4.45) 



where we have used the fact that in coordinates along j^{(f> o (fix 1 ), 



Integrating ( 4.45 ) by parts, we obtain 



u x 



91 (j^o^-^^-Ci^o^)) 



dy AKJ VY " ^ " dx^ 



V A d n+1 x 



+ 



dL 

dv A „ 



^"^l ^I^V fiJ. (4.46) 



<9x 



a 



Let a be the n-form integrand of the boundary integral in ( 4.46| ); then Jg Ux 
fdj 1 (<j>o<j>~ 1 )(u x ) s i nce a * s invariant under this lift. Additionally, from equation ( |4.39| j, 
we obtain the horizontal contribution 



dSj, ■ V h = [ (V^) J (Ld n+1 x) = [ 

JdU x Jd, 



V^Ld n X 



dj^^X-^Ux) 



I" 



(4.47) 



so combining equations ( 4.46j ) and ( [4.47 ), a simply computation verifies that 

d dL 



dSs-V 



Ux 



dT 



Q x fi Q V A 



-(/(^o^- 1 )) 



d n+1 x ®dy A \ - V 



+ 



V 



( Pit 

^{^A-U'^^X^Ad-X, 



+ 



L 



dL 

dv A , 



(J )) «ZI 



xY 



dxt 1 



d n+i x } . (4.48) 



The vector V in the second term of ( 4,4§| ) may be replaced by j 1 (V) since t^y,j 1 (y)- 
vertical vectors are clearly in the kernel of the form that V is acting on. This shows 
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that (4.43) and ( |4.44| ) hold, and hence that the boundary integral in ( |4.48| ) may be 



written as 



IdUx 

Now, if we choose another coordinate chart O2, the coordinate expressions of 



Del£ and 0£ must agree on the overlap 0± n O2 since the left-hand-side of ( |4.4l|) 
is intrinsically defined. Thus, we have uniquely defined Del£ and 0£ for any Ux 
such that U x PidX = $. 



Finally, ( [4.42 ) holds, since Qc = d&c is also intrinsically defined and both 



sides of the equation yield the same coordinate representation, the Euler-Lagrange 
equations in Ux- B 



Remark To prove Theorem |4,4| for the case Ux = X, we must modify the proof 
to take into account the boundary conditions which are prescribed on dX. 

Corollary 4.2 The (n + l)-form ®c defined by the variational principle satisfies 
the relationship 

C(j 1 (z)) = z*@c 
for all holonomic sections z G T('^x,j 1 (y))- 



Proof. This follows immediately by substituting (4.42) into ( |4.41| ) and integrating 
by parts using Cartan's formula. ■ 

Remark We have thus far focused on holonomic sections of J^Y), those that are 
the first jets of sections of Y, and correspondingly, we have restricted the general 
splitting of TY given by 

TY = image 7 © VY for any 7 G T(J 1 (y)), 



to TY = Tcp © VY, (ft G r(Y) as we specified in Q4.38 ). For general sections 



7 G r(J 1 (y)), the horizontal bundle is given by image 7, and the Frobenius theorem 
guarantees that 7 is locally holonomic if the connection is flat, or equivalently if the 
curvature of the connection i? 7 vanishes. Since this is a local statement, we may 
assume that Y = U x R , where U C M n+1 is open, and that ttxy is simply the 
projection onto the first factor. For <fi G T(Y), and 7 G F(J 1 (Y)), 7(2:, <f>( x )) '■ 
jjn+i j^v j g a ii near operator which is holonomic if <p'(x) = j(x, (j>(x)), where 
<j)'{x) is the differential of (ft, and this is the case whenever the operator <ft"(x) is 
symmetric. Equivalently, the operator 

S y (x, y) ■ (v, w) := £>i7(z, y) ■ {v, w) + ^27(2;, V) ■ {l{x, y) ■ v, w) 

is symmetric for all v,w G M n+1 . One may easily verify that the local curvature is 
given by 

Rj{x,y) ■ (v,w) := S^(x,y) ■ (v,w) - 5 7 (x,y) • (w,v) 
and that 7 = j l (cft) locally for some (ft G T(Y), if and only if i? 7 = 0. 
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The variational route to Noether's Theorem. Suppose the Lie group Q acts 
on C and leaves the action S invariant so that 



«S($(7/y») = S(</>) for all rjy 6 Q. (4.49) 

This implies that for each rjy € G, $>(j]Y,<j>) G T 3 whenever eft € V. We restrict 
the action of C/ to T 3 , and let £c be the corresponding infinitesimal generator on C 
restricted to points in V; then 

o = (e c jd5) = / /OWxW^je^] 



since %i(#)0£ = by (pUSD and Corollary [D 



We denote the covariant momentum map on J 1 (1 A ) by 6 L(g, A"(J 1 (y)) 
which we define as 

j l {O^Sl c = dJ c {t). (4.50) 

Using ( p0| ), we find that f D d\j l {4> o ^)*J C {C)] = 0, and since this must 
hold for all infinitesimal generators £e a f </• G C, the integrand must also vanish so 
that 

d[jV°0xr^(O]=O, (4.51) 
which is precisely a restatement of the covariant Noether Theorem. 



5 Veselov-type Discretizations of Multisymplectic Field 
Theory 

5.1 General Theory 

We now generalize the Veselov discretization given in Section (3) to multisymplectic 
field theory, by discretizing the spacetime X. For simplicity we restrict to the 
discrete analogue of dimX = 2; i.e. n = 1. Thus, we take I = ZxZ = {(hj)} an d 
the fiber bundle Y to be X x T for some smooth manifold T . 



Notation. The development in this section is aided by a small amount of notation 
and terminology. Elements of Y over the base point are written as and the 
projection ttxy acts on Y by i^xYiyij) = (hj)- The fiber over E X is denoted 
Yij. A triangle A of X is an ordered triple of the form 

A= + + + 1)). 

The first component of A is the first vertex of the triangle, denoted A 1 , 

and similarly for the second and third vertices. The set of all triangles in X is 
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denoted X . By abuse of notation the same symbol is used for a triangle and the 
(unordered) set of its vertices. A point £ X is touched by a triangle if it is 
a vertex of that triangle. If U C X, then € U is an interior point of U if 
U contains all three triangles of X that touch The interior intU of U is 

the collection of the interior points of U. The closure clU of U is the union of all 
triangles touching interior points of U. A boundary point of U is a point in U 
and cl U which is not an interior point. The boundary of U is the set of boundary 
points of U, so that 

dU=(Un cl U) \ int U 

Generally, U properly contains the union of its interior and boundary, and we call 
U regular if it is exactly that union. A section of Y is a map <f> : U C X — > 1" 
such that 7rxy o = idjy . 




Figure 5.1: Depiction of the heuristic interpretation of an element of J X Y when X is 
discrete. 

Multisymplectic phase space. We define the first jet bundl^ of Y to be 

jly = {{yij,Vij+i,yi+ij+i) | e X, yij,yij+i,yi+ij+i e 
= X A x^ 3 . 

Heuristically (see Figure ( |5.1| )), A corresponds to some grid of elements xij in con- 
tinuous spacetime, say X, and (yij, Vij+i, y^ij+x) G J 1 !^ corresponds to j l <f){x), 
where x is "inside" the triangle bounded by Xij,Xij+i,Xi+ij+x, and <j> is some smooth 
section of X x T interpolating the field values yij,yij+i, yi+i j+i- The first jet ex- 
tension of a section of y is the map j l cf):X A ^J l Y defined by 

jV(A) = (A, ^A 1 ), 0(A 2 ), <p(A 3 )) . 

Using three vertices is the simplest choice for approximating the two partial derivatives of the 
field (j>, but may not lead to a good numerical scheme. Later, we shall also use four vertices together 
with averaging to define the partial derivatives of the fields. 
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Given a vector field Z on Y, we denote its restriction to the fiber Y\j by Zjj, and 
similarly for vector fields on J X Y . The first jet extension of a vector field Z on 
V is the vector field j l Z on J l Y defined by 

j 1 Z{y A i,y A 2,y A s) = (Z A i(y A i), Z A 2(y A i), Z A3 (y A3 )) , 

for any triangle A. 







(i + hj) 


0+1,7+ 1) 




A 


A 




(U-i) 


A 


a, j) 


(U+i) 


- 1,7- 1) 




(i-hj) 





Figure 5.2: TTie triangles which touch 



The variational principle. Let us posit a discrete Lagrangian L : J l Y — > R. 

Given a triangle A, define the function L A : J^" 3 — > R by 

LA(yi,y2,y3) = L(A,y 1 ,y 2 ,y3), 

so that we may view the Lagrangian L as being a choice of a function L A for each 
triangle A of X. The variables on the domain of L A will be labeled y 1 ,?/ 2 ,?/ 3 , 
irrespective of the particular A. Let E/ be regular and let Cjj be the set of sections 
of Y on U, so Cu is the manifold J^I^L The action will assign real numbers to 
sections in Cu by the rule 

5(0)= ^ LojV(A). (5.1) 

A;AC{/ 

Given <fi £ Cu and a vector field V, there is the 1-parameter family of sections 

(FYc(>)(i,j)=FY ij mj)), 

where F Vi ' denotes the flow of Vij on T. The variational principle is to seek 
those 6 for which 



d_ 
~de 



S(FU) = 

e=0 



for all vector fields V. 
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The discrete Euler-Lagrange equations. The variational principle gives cer- 
tain field equations, the discrete Euler-Lagrange field equations (DELF equa- 
tions), as follows. Focus upon some (i,j) G intC/, and abuse notation by writing 
(f>(i,j) = yij. The action, written with its summands containing y^ explicitly, is (see 
Figure Q) 

S = h L(yij,yij + i,y i+ ij + i) + L(y ij - 1 , yy, yi+ij) + L(yt-ii-i,yi-ii,2/tj) H 

so by differentiating in y^, the DELF equations are 
dL . dL dL 

^-r(yii,2/ii+i,yi+ii+i) + ^(i/ii-i^w+ii) + ^3 (i/i-ii-i.w-i^i/ii) = 0, 

for all (i,j) € intf/. Equivalently, these equations may be written 

J;A;(t,i)=A' 

for all (£, j) G int C7. 

The discrete Cartan form. Now suppose we allow nonzero variations on the 
boundary dU, so we consider the effect on S of a vector field V which does not 
necessarily vanish on dU. For each G 9C/ find the triangles in U touching 

There is at least one such triangle since G clU; there are not three such 
triangles since (i,j) intt/. For each such triangle A, occurs as the Z th vertex, 
for one or two of I = 1,2,3, and those Z th expressions from the list 

Jjjjr (yij ,yij+i,yi+ij+i) Vij [yij ) , 
if^iyt j-u Vij,yi+i j)Vij(yij), 

dL 
dy 3 

yielding one or two numbers. The contribution to dS from the boundary is the 
sum of all such numbers. To bring this into a recognizable format, we take our cue 
from discrete Lagrangian mechanics, which featured two 1-forms. Here the above 
list suggests the three 1-forms on J l Y , the first of which we define to be 

Q l L (y lj ,y ij+1 ,y i+ i j+1 ) • (%^% J+1 ,% +lj+1 ) 



^-3 (yi-ij-i,yi-ij,yij) Vij (y^ ) 



dL 

= Q-T{yij>yij+i>yi+ij+i) • (% j5 o,o), 

6^ and 0| being defined analogously. With these notations, the contribution to dS 
from the boundary can be written 6l(^>) ■ V, where 6l is the 1-form on the space of 
sections Cjj defined by 

em-v= E I £ [(iVr(/^jeL)](A)] . (5.3) 

A;An9C/^0 \l;A l edU ' J 
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In comparing (^^) with ( |4.41| ), the analogy with the multisymplectic formalism of 
Section (4) is immediate. 



The discrete multisymplectic form formula. Given a triangle A in X, we 
define the projection 7r A : Cy — > J l Y by 

tta(0) = (A,y A i,y A 2,y A 3). 

In this notation, it is easily verified that (^^) takes the convenient form 

0l= Yl f E ^©U- (5-4) 

A;An<9£/^0 \l;A l £dU J 

A first-variation at a solution ^ of the DELF equations is a vertical vector field V 
such that the associated flow F v maps 4> to other solutions of the DELF equations. 
Set Q l L = -dQ l L . Since 

@l + e| + e| = dL, (5.5) 

one obtains 

nl + n 2 L + n| = o, 

so that only two of the three 2-forms Q l L , I = 1, 2, 3 are essentially distinct. Exactly 
as in Section (2), the equation d 2 S = 0, when specialized to two first-variations V 



and W now gives, by taking one exterior derivative of (5.4) 



Q = d0 L {4>){V,W)= E VjWJttWl 

A;AnO(7^0 \l;A l &dU 



which in turn is equivalent to 



E ( E [tfwtfvjfwjt^w] =o. 

A;And£/^0 \l;A l edU J 



(5.6) 



Again, the analogy with the multisymplectic form formula for continuous space- 
time ( [4.18| ) is immediate. 



The discrete Noether theorem. Suppose that a Lie group G with Lie algera g 
acts on F by vertical symmetries in such a way that the Lagrangian L is G- invariant. 
Then G acts on Y and J l Y in the obvious ways. Since there are three Lagrange 1- 
forms, there are three momentum maps J\ I = 1,2,3, each one a g*-valued function 
on triangles in X, and defined by 

4 = Oiy- 1 ©l> 
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for any £ E g. Invar iance of L and 



imply that 



J 1 + J 2 + J 3 



0. 



so, as in the case of the 1-forms, only two of the three momenta are essentially 
distinct. For any £, the infinitesimal generator £y is a first-variation, so invariance 
of S, namely £y J dS = , becomes £y J 0l = 0. By left insertion into (5J5), this 
becomes the discrete version of Noether's theorem: 



E E Aa) 



(5.7) 



A;An<9t/^0 \l;A l £dU 




Time 



Space 



Figure 5.3: Symplectic flow and conservation of momentum from the discrete Noether 
theorem when the spatial boundary is empty and the temporal boundaries agree. 



Conservation in a space and time split. To understand the significance of ( p.6|) 
and ( |5.7| ) consider a discrete field theory with space a discrete version of the circle 
and time the real line, as depicted in Figure (|5,3j), where space is split into space 
and time, with "constant time" being constant j and the "space index" 1 < i < N 
being cyclic. Applying (5/7) to the region \ j = 0,1,2} shown in the Figure, 

Noether's theorem takes the conservation form 



N 



N 



E jl (yio,ya,yi+ii) = -E( j2 ( yil > yi2,2/i + 12 ) + j3 (y*i' ^2,^+12)) 



8=1 



i=l 



N 



i=l 

Similarly, the discrete multisymplectic form formula also takes a conservation form. 
When there is spatial boundary, the discrete Noether theorem and the discrete 
multisymplectic form formulas automatically account for it, and thus form nontrivial 
generalizations of these conservation results. 
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Furthermore, as in the continuous case, we can achieve "evolution type" sym- 
plectic systems (i.e. discrete Moser-Veselov mechanical systems) if we define Q as 
the space of fields at constant j, so Q = J- , and take as the discrete Lagrangian 

N 
i=l 

Then the Moser-Veselov DEL evolution-type equations (|3.2| ) are equivalent to the 
DELF equations (|5.2| ), the multisymplectic form formula implies symplecticity of 
the Moser-Veselov evolution map, and conservation of momentum gives identical 
results in both the "field" and "evolution" pictures. 



Example: nonlinear wave equation. To illustrate the discretization method 
we have developed, let us consider the Lagrangian ( |Og ) of Section (4), which 
describes the nonlinear sine-Gordon wave equation. This is a completely integrable 
system with an extremely interesting hierarchy of soliton solutions, which we shall 
investigate by developing for it a variational multisymplectic-momentum integrator; 
see the recent article by Palais [1997] for a wonderful discussion on soliton theory. 

To discretize the continuous Lagrangian, we visualize each triangle A as having 
base length h and height k, and we think of the discrete jet (i/a 1 ; 2/A 2 > 2/A 3 ) as 
corresponding to the continuous jet 

d4> ,_ , _ Vij+i - Vjj dcf> _ y i+ i j+1 - yi j+1 
dx* {Vij) ~ h ' dx l[Vi3 '~ k 

where t/y is a the center of the triangle This leads to the discrete Lagrangian 

L = i ( V2-yi \ 2 l / m -V2 \ 2 N ( yi + y2 + y3 

2\h) 2\k) 
with corresponding DELF equations 

Vi+ij -2yij + Vi-ij _ Vij+i -tyij + yij-i 

k 2 h 2 



, 1 N > ( Vis ±yi±j± + yi+w \ 

3 v 3 ; 

, i jy/ ( Vij-i + yij + yi+ij \ 

3 v 3 ; 



+ In' ( y^+y^i+y^ = o. ( 5. 8 ) 

When N = (wave equation) this gives the explicit method 

h 2 , 

y* J'+l = If 3 ~ 2 Vi3 + 3) + 2 Vij ~ J-l' 

which is stable whenever the Courant stability condition is satisfied. 

2 Other discretizations based on triangles are possible; for example, one could use the value j/y 
for insertion into the nonlinear term instead of yij. 
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Extensions: Jets from rectangles and other polygons. Our choice of discrete 
jet bundle is obviously not restricted to triangles, and can be extended to rectangles 
or more general polygons (left of Figure( |5.4| )). A rectangle is a quadruple of the 
form, 

A = ((*, j), + 1), (« + 1,3 + 1), (i + l,i)), 

a point is an interior point of a subset U of rectangles if U contains all four rect- 
angles touching that point, the discrete Lagrangian depends on variables yi, • • • ,2/4, 
and the DELF equations become 

dL dL 

i^(yijiyij+i,yi+ij+i,yi+ij) + -Q-2\yij-i,yij,yi+ij,yi+ij-i) 

dL . . dL . 

+Q-3{yi-ij-i,yi-ij,yij,yij-i) + ^w-i ^-1^+1,^+1,%) = o. 

The extension to polygons with even higher numbers of sides is straightforward; one 
example is illustrated on the right of Figure( |5.4|). The motivation for consideration 



(i+lj-l) 




a +1,7) 




UJ-i) 


i 

v , 
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v 










(i-lj-l) 












(i-hj) 








Figure 5.4: On i/ie Ze/i, £/ie method based on rectangles; on the right, a possible method 
based on hexagons. 

of these extensions is enhancing the stability of the triangle-based method in the 
nonlinear wave example just above. 



Example: nonlinear wave equation, rectangles. Think of each rectangle A 
as having length h and height k, and each discrete jet {y^-,y&2-,y/\3-,y/\±) being 
associated to the continuous jet 

d(p , , Vij+i ~ y%j dcj) , . If y i+lj -y i:j Vi+ij+x - Vij+i\ 
W {P) = h ' ^ P) = 2{ k + k J' 

where p is a the center of the rectangle. This leads to the discrete Lagrangian 

2 V h J 2\ 2k 2k J 
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N 



V1+V2 + V3+ 2/4 



(5.9) 



If, for brevity, we set 



Vij + Vi 3+1 + Vi+l j+1 + Vi+li 



then one verifies that the DELF equations become 



1 y i+ ij -2yjj + yj-ij 1 Vi+i j+i - 2y% j+i + Vi-i j+i 



k 2 



+ 



A: 2 



4 A; 2 



gfcj+i - ^ Vij + yij-i 
/i 2 



+ 



N'iyij) + iV'fe^i) + N'fa..!^) + iV'G/i-^) 



which, if we make the definitions 

V2 



o> fe ytj = Vi j+i ~ 2 Vij + Vi j-i, d k yij = y i+1 3 -2 y i3 - + j/j_i j , 



is (more compactly) 



1 



~d%y ij+ i + ^dlyij + jdlyij-x 



j^d 2 hyij + TV'(^) = 0. 



(5.10) 



These are implicit equations which must be solved for yij+i, 1 < % < iV, given yjj, 
Vij-ii 1 < * < -W; rearranging, an iterative form equivalent to ( 5.10| ) is 



/i 2 



2(/i 2 + 2A; 2 ) 



h 2 



2{h 2 + 2k 2 ) 



Ui-lj+l 



h 2 



h 2 + 2A; 2 



([Vi+ij ~ 2 Vij +Vi-lj) 



+ 
+ 



2A; 2 



h 2 + 2k 2 
h 2 k 2 



2{h 2 + 2k 2 ) 



+ -^{Vi+ij-l -^Vij-i + Vi-lj-i] 



{2 yij - yij-i) 

(N'iyij) + N'faj-i) + N'im-ij-x) + N'iy^j)). 



In the case of the sine-Gordon equation the values of the field ought to be considered 
as lying in § , by virtue of the vertical symmetry y ^ y + 2-k. Soliton solutions 
for example will have a jump of 2tt and the method will fail unless field values at 
close-together spacetime points are differenced modulo 2ir. As a result it becomes 
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important to calculate using integral multiples of small field-dependent quantities, 
so that it is clear when to discard multiples of 2tt, and for this the above iterative 
form is inconvenient. But if we define 

dhVij = Vij+1 ~ Viji dkDij = Ui+lj ~ Viji 

then there is the following iterative form, again equivalent to ( |5.10| ) 

Vij+i = Vij + dlyij, and 



h 2 



2(h 2 + 2k 2 ) 
h 2 



h 2 + 2k 2 

h 2 k 2 

+ 



2k 2 



i-lj 



2{h 2 + 2k 2 ) 



(5.11) 



One can also modify (^)1]) so as to treat space and time symmetrically, which 
leads to the discrete Lagrangian 



L 



2 V 2h 2h 



1 ( V4 - m V3 - V2 

2 \ 2k " 2k 



+N 



V1+V2+V3 + 2/4 



and one verifies that the DELF equations become 



1 

fc2 



1 

h 2 



N'(y, 



i.i' 



an equivalent iterative form of which is 

and 

dhVi+ij + dhVij 



Vij+i = Vij + d\y, 
h 2 -k 2 



2{h 2 + k 2 ) 
h 2 



h 2 -k 2 
2(h 2 + k 2 ) 



dhVi-lj 



2{h 2 + k 2 ) 
h 2 



(3d 2 y l3 + tfvij-!) 



+ 



2{h 2 + k 2 ) 
h 2 k 2 



( 2d hVij + d hVi+ij +dhVi- 



+ 



2{h 2 + k 2 ) 



(5.12) 



(5.13) 
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Figure 5.5: Top left: the wave forms for a two soliton kink and antikink collision us- 
ing (5.1i). Top right: the energy error. Bottom left: the wave form at time t sa 11855. 
Bottom right: the portion of the bottom left graph for spatial grid points 1 ... 16. 



5.2 Numerical checks. 

While the focus of this article is not the numerical implementation of the integrators 
which we have derived, we have, nevertheless, undertaken some preliminary numeri- 
cal investigations of our multisymplectic methods in the context of the sine-Gordon 
equation with periodic boundary conditions. 



The rectangle-based multisymplectic method. The top half of Figure ( |5,5D 
shows a simulation of the collision of "kink" and "antikink" solitons for the sine- 
Gordon equation, using the rectangle-based multisymplectic method ( 5.12| ). In the 
bottom half of that figure we show the result of running that simulation until the 
solitons have undergone about 460 collisions; shortly after this the simulation stops 
because the iteration ( 5.13] ) diverges. The anomalous spatial variations in the wave- 
form of the bottom left of Figure ( |5.5[) have period 2 spatial grid divisions and 
are shown in finer scale on the bottom right of that figure. These variations are 
reminiscent of those found in Ablowitz, Herbst and Schober [1996] for the com- 
pletely integrable discretization of Hirota, where the variations are attributed to 
independent evolution of waveforms supported on even vs. odd grid points. Obser- 
vation of ( |5.12 ) indicates what is wrong: the nonlinear term N contributes to ( 5.12 ) 
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in a way that will average out these variations, and consequently, once they have 
begun, ( |5.12| ) tends to continue such variations via the linear wave equation. In 
Ablowitz et. al., the situation is rectified when the number of spatial grid points 
is not even, and this is the case for ( |5.12| ) as well. This is indicated on the left 
of Figure ( |5.6| ), which shows the waveform after about 5000 soliton collisions when 
N = 255 rather than N = 256. Figure ( |5.7| ) summarizes the evolution of energy 
error^j for that simulation. 

Initial data. For the two-soliton-collision simulations, we used the following initial 
data: h = k/8 (except h = k/16 where noted), where k = 40/iV and N = 255 spatial 
grid points (except Figure (|5.5[ ) where N = 256). The circle that is space should 
be visualized as having circumference L = 40. Let k = 1 — e where e = 10 -6 , 
L = L/4 = 10, 



and 



/•!/« 1 / L 2 

P = 2 _ dy « 15.90, c = \ 1 - T ra .7773, 



2 arcsin ( sn ( , ; k 



c 



2 



Then (f>(x — ct) is a kink solution if space has a circumference of L. This kink and 
an oppositely moving antikink (but placed on the last quarter of space) made up 
the initial field, so that t/io = </>(40(i — 1)/N), i = 1, . . . ,N, where 

(<j>{x) < x < L/4 

2vr L/4 < x < 3L/4 , 

2vr-</>(x-3L/4) 3L/4 < x < L 

while yn = yio + 0(4O(z — 1)/N)h where 

(4>(x - he) - <p{x))/h < x < L/4 
cj)(x) = { L/4 < x < 3L/4 . 

-((f){x - he) - c/){x))/h 3L/4<x<L 



The discrete energy that we calculated was 

JV 



E 



IfVij+l—yij , Vi+lj+l — Vi+1 j 



2 V 2h 2h 
+ 



2 



1 f yi+ij -Vij , yi+ij+i - yij+i \ 2 ATl - A 

2 I 2k + 2k ) ly ^>J- 
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Figure 5.6: On the left, the final wave form (after about 5000 soliton collisions at t w 
129133,) obtained using the rectangle-based multisymplectic method (5.1i). On the right, 
the final waveform (at t rs 129145 ^ from the energy- conserving method (5.14) of Vu-Quog 
and Li. In both simulations, temporal drift is occurring. For this reason the waveforms are 
inverted with respect to one another; moreover, the separate solitons are drifting at slightly 
different rates, as indicated by the off-center waveforms. 
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Figure 5.7: On the left, the energy error corresponding to our multisymplectic method 



(5. IS) for 5000 solition collisions; the three graphs correspond to the minimum, average, 
and maximum energy error over consecutive 5000 time step regions. On the right, the final 
energy error (i.e. the energy error after about 5000 soliton collisions), which can be compared 
with the initial energy error plot in the top left of Figure (pT, 
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Comparison with energy-conserving methods. As an example of how our 
method compares with an existing method, we considered the energy-conserving 
method of Vu-Quoc and Li [1993], page 354: 



1 



h 2 



dhij 



+ 



1 f N(y ij+1 )-N( yij ) | Njy^-Njy^^) 

2 V Vi j+l ~ Vij Vij ~ Vi j-i 



0. (5.14) 



This has an iterative form similar to ( 5.13 ) and is quite comparable with ( 5.10 ) 
and ( |5,12| ) in terms of the computation required. Our method seems to preserve the 
soliton waveform better than ( 5.14 ), as is indicated by comparison of the left and 
right Figure (|5.6|) . 

In regards to the closely related papers Vu-Quoc and Li [1993] and Li and Vu- 
Quoc [1995], we could not verify in our simulations that their method conserves 
energy, nor could we verify their proof that their method conserves energy. So, as 
a further check, we implemented the following energy-conserving method of Guo, 
Pascual, Rodriguez, and Vazquez [1986]: 



d kVij 



\Vij + 

which conserves the discrete energy 

yij)(yij — Vij- 



N( yi 



N(y 



Vij+1 — Vij-1 



(5.15) 




h? 



±> + i 

2 



Vi+i j ~ Vij 
k 



N(y. 



This method diverged after just 345 soliton collisions. As can be seen from ( 5.15| ), 
the nonlinear potential ./V enters as a difference over two grid spacings, which sug- 
gests that halving the time step might result in a more fair comparison with the 
methods ( 5.12j ) or ( |5.14| ), With this advantage, method ( |5.15| ) was able to simulate 
5000 soliton collisions, with a waveform degradation similar to the energy-conserving 
method ( 5.14| ), as shown at the bottom right of Figure ( |5.8| ). The same figure also 



shows that, although the energy behavior of (5.15) is excellent for short time simula- 
tions, it drifts significantly over long times, and the final energy error has a peculiar 
appearance. Figure ( |5.9D shows the time evolution of the waveform through the 
soliton collision that occurs just before the simulation stops. Apparently, at the 
soliton collisions, significant high frequency oscillations are present, and these are 
causing the jumps in the energy error in the bottom left plot of Figure (|5.8| ). This 
error then accumulates due to the energy-conserving property of the method. In 
these simulations, so as to guard against the possibility that this behavior of the 



energy was due to inadequately solving the implicit equation ( 5.15 ), we imposed a 
minimum limit of 3 iterations in the corresponding iterative loop, whereas this loop 
would otherwise have converged after just 1 iteration. 
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Figure 5.8: A long-time simulation using the energy- conserving method (5.11 ) of Guo et al. 
Above left: the initial energy error. Above right: the average energy error over consecutive 
5000 time step regions (the maximum and minimum closely parallel the average). Below 
left: the final energy error. Below right: the final waveform at t « 129149. 



Comparison with the triangle-based multisymplectic method. The dis- 



crete second derivatives in the method (5.15) are the same as in the triangle- 



based multisymplectic method (|5.8j ); these derivatives are simpler than either our 
rectangle-based multisymplectic method ( |5.12D or the energy-conserving method of 
Vu-Quoc and Li ( 5.14| ). To explore this we implemented the triangle-based multi- 



symplectic method ( |5.8| ) . Even with the less complicated discrete second derivatives 
our triangle-based multisymplectic method simulated 5000 soliton collisions with 
comparable energy ^ and waveform preservation properties as the rectangle-based 
multisymplectic method ( |5,12| ), as shown in Figure fl5,ll| ). Figure ( |5.10| ) shows the 
time evolution of the waveform through the soliton collision just before the sim- 
ulation stops, and may be compared to Figure ( |5.9| ). As can be seen, the high 
frequency oscillations that are present during the soliton collisions are smaller and 
more smooth for the triangle-based multisymplectic method than for the energy- 

4 The discrete energy that we calculated was 



N / 1 

i=l v 
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Figure 5.9: T/ie soliton collision at time t w 129130, a/ter £/ie energy- conserving 
method (5.11) of Guo et al. has simulated about 5000 soliton collisions. The solitons collide 



beginning at the top left and proceed to the top right, then to the bottom left, and finally 
to the bottom right. The vertical scales are not constant and visually exaggerate the high 
frequency oscillations, which are small on the scale to 2ir. 
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Figure 5.10: Similar to the above plot but for our triangle-based multisymplectic 
method (5.1 ). 
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conserving method ( 5.15j ). A similar statement is true irrespective which of the two 
multisymplectic or two energy conserving methods we tested, and is true all along 
the waveform, irrespective of whether or not a soliton collision is occurring. 
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Figure 5.11: A simulation of 5000 soliton collisions using the triangle-based multisym- 
plectic method (5.t ). Above left: The initial energy error. Above right: The minimum, 
average and maximum energy as in the left of Figure \5.\ ). Below left: the final waveform 
(at t rs 129130,). Below right: the final energy error. 



Summary. Our multisymplectic methods are finite difference methods that are 
computationally competitive with existing finite difference methods. Our methods 
show promise for long-time simulations of conservative partial differential equations, 
in that, for long-time simulations of the sine-Gordon equation, our method 1) had 
superior energy-conserving behavior, even when compared with energy- conserving 
methods] 2) better preserved the waveform than energy-conserving methods; and 3) 
exhibited superior stability, in that our methods excited smaller and more smooth 
high frequency oscillations than energy-conserving methods. However, further nu- 
merical investigation is certainly necessary to make any lasting conclusions about 
the long-time behavior of our integrator. 

The programs. The programs that were used in the preceding simulations are 
"C" language implementations of the various methods. A simple tridiagonal LUD 
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method was used to solve the linear equations (e.g. the left side of ( 5.13; )), as in 
Vu-Quoc and Li [1993], page 379. An 8 th order extrapolator was used to provide 
a seed for the implicit step. All calculations were performed in double precision 
while the implicit step was terminated when the fields ceased to change to single 
precision; the program's output was in single precision. The extrapolation usually 
provided a seed accurate enough so that the methods became practically explicit, in 
that for many of the time-steps the first or second run through the iterative loops 
solving the implicit equations solved those equations to single precision. However, 
in the absence of a regular spacetime grid the expenses of the extrapolation and 
solving the linear equation would grow. Our programs are freely available at URL 
http : //www. cds . caltech.edu/cds. 



6 Concluding Remarks 

Here we make a few miscellaneous comments and remark on some work planned for 
the future. 



Lagrangian reduction. As mentioned in the text, it is useful to have a discrete 
counterpart to the Lagrangian reduction of Marsden and Scheurle [1993a, b], Holm, 
Marsden and Ratiu [1998a] and Cendra, Marsden and Ratiu [1998]. We sketch 
briefly how this theory might proceed. This reduction can be done for both the case 
of "particle mechanics" and for field theory. 

For particle mechanics, the simplest case to start with is an invariant (say left) 
Lagrangian on the tangent bundle of a Lie group: L : TG — > R. The reduced 
Lagrangian is I : g — > R and the corresponding Euler-Poincare equations have a 
variational principle of Lagrange d'Alembert type in that there are constraints on 
the allowed variations. This situation is described in Marsden and Ratiu [1994]. 

The discrete analogue of this would be to replace a discrete Lagrangian L : 
G x G — > R by a reduced discrete Lagrangian £ : G — * R related to L by 

Kai92 l ) = L (#i>52) 

In this situation, the algorithm from G x G to G x G reduces to one from G to 
G and it is generated by £ in a way that is similar to that for L. In addition, the 
discrete variational principle for L which states that one should find critical points 
of 

with respect to gi to implicitly define the map (51,52) l— ► (92,93), reduces naturally 
to the following principle: Find critical points of 1(g) +l(h) with respect to variations 
of g and h of the form g£ := L 5 £ and £h = Rh£ where L g and Rh denote left and 
right translation and where £ € 0. In other words, one sets to zero, the derivative of 
the sum £(gg~ 1 ) + i(g e h) with respect to e at e = for a curve g e in G that passes 
through the identity at e = 0. This defines (with caveats of regularity as before) 
a map of G to itself, which is the reduced algorithm. This algorithm can then 
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be used to advance points in G x G itself, by advancing each component equally, 
reproducing the algorithm on G x G. In addition, this can be used with the adjoint 
or coadjoint action to advance points in q or q* to approximate the Euler-Poincare 
or Lie-Poisson dynamics. 

These equations for a discrete map, say (f>e : G — > G generated by I on G are called 
the discrete Euler Poincare equations as they are the discrete analogue of the 
Euler-Poincare equations on q. Notice that, at least in theory, computation can be 
done for this map first and then the dynamics on G x G is easily reconstructed by 
simply advancing each pair as follows: (51,52) l— ► (%i,%2), where h = 4>i(g\ <fe)- 

If one identifies the discrete Lagrangians with generating functions (as explained 
in Wendlandt and Marsden [1997]) then the reduced Lagrangian generates the re- 
duced algorithm in the sense of Ge and Marsden [1988] , and this in turn is closely 
related to the Lie-Poisson-Hamilton-Jacobi theory. 

Next, consider the more general case of TQ with its discretization Q x Q with 
a group action (assumed to be free and proper) by a Lie group G. The reduction 
of TQ by the action of G is TQ/G, which is a bundle over T(Q/G) with fiber 
isomorphic to g. The discrete analogue of this is (Q x Q)/G which is a bundle 
over (Q/G) x {Q/G) with fiber isomorphic to G itself. The projection map tt : 
(Q x Q)/G -»■ (Q/G) x (Q/G) is given by [(91,92)] ^ ([Qi], [92]) where [ ] denotes 
the relevant equivalence class. Notice that in the case in which Q = G this bundle 
is "all fiber" . The reduced discrete Euler-Lagrange equations are similar to those in 
the continuous case, in which one has shape equations couples with a version of the 
discrete Euler-Poincare equations. 

Of course all of the machinery in the continuous case can be contemplated here 
too, such as stability theory, geometric phases, etc. In addition, it would be useful 
to generalize this Lagrangian reduction theory to the multisymplectic case. All of 
these topics are planned for other papers. 



Role of uniformity of the grid. Consider an autonomous, continuous La- 
grangian C : TQ — ► M where, for simplicity, Q is an open submanifold of Euclidean 
space. Imagine some not necessarily uniform temporal grid (to, t\, • • • ) of R, so that 
to < t\ < t2 < ■ ■ ■ ■ In this situation, it is natural to consider the discrete action 

fe=i fc=i ^ k k-\j 

This action principle deviates from the action principle (|3.l| ) of Section 3 in that the 
discrete Lagrangian density depends explicitly on k. Of course nonautonomous con- 
tinuous Lagrangians also yield ^-dependent discrete Lagrangian densities, irrespec- 
tive of uniformity of the grid. Thus, nonuniform temporal grids or nonautonomous 
Lagrangians give rise to discrete Lagrangian densities which are more general those 
those we have considered in Section (3). For field theories, the Lagrangian in the 
action fl5.1| ) depends on the spacetime variables already, through its explicit depen- 
dence on the triangle A. However, it is only in the context of a uniform grid that we 
have experimented numerically and only in that context that we have discussed the 
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significance of the discrete multisymplectic form formula and the discrete Noether 
theorem. 

Using ([O]) as an example, will now indicate why the issue of grid uniformity 
may not be serious. The DEL equations corresponding to the action ( |6.l] ) are 

^(gfc, gfc-i) + ^o fc+1 (gfc+l, gfc) = 0, k=l,2, (6.2) 
oq\ dq 2 

and this gives evolution maps Ffc+i.fc : Q X Q — ► Q X Q defined so that 

Fk+i,k(qk,qk-i) = (qk+i,qk), k = 1,2, • • • 

when ( |6.2[) holds. For the canonical 1-forms corresponding to (|3.4j ) and ( |3.5[ ) we 
have the /^-dependent one forms 

dL 

OZ k (qi,<lo) ■ (Sqi,Sq ) = -^—(qi,qo)$qo, (6.3) 



and 



dqo 



dL 

?L,fc(9i>?o) • (<^i,<fyo) = -3— (?i,9o)*?i 3 (6-4) 



and Equations ( |3.7D and ( |3.9| ) become 

Fk+x,k( dd tk) = -MZj+v de I,k + de U = ® (6-5) 
respectively. Together, these two equations give 

^ + i,k( d0 ik) = <,k+v (6-6) 

and if we set 

Fk = Fk,k-l Fk-i t k-2 o • • • o F 2j i 

then (|6.6[) chain together to imply F^dO^) = d0^ k . This appears less than 
adequate since it merely says that the pull back by the evolution of a certain 2-form 
is, in general, a different 2-form. The significant point to note, however, is that this 
situation may be repaired at any k simply by choosing = L\. It is easily verified 
that the analogous statement is true with respect to momentum preservation via 
the discrete Noether theorem. 

Specifically, imagine integrating a symmetric autonomous mechanical system in 
a timestep adaptive way with Equations fl6.2|) . As the integration proceeds, various 
timesteps are chosen, and if momentum is monitored it will show a dependence on 
those choices. A momentum-preserving symplectic simulation may be obtained by 
simply choosing the last timestep to be of equal duration to the first. This is the 
highly desirable situation which gives us some confidence that grid uniformity is 
a nonissue. There is one caveat: symplectic integration algorithms are evolutions 
which are high frequency perturbations of the actual system, the frequency being 
the inverse of the timestep, which is generally far smaller than the time scale of any 
process in the simulation. However, timestep adaptation schemes will make choices 
on a much larger time scale than the timestep itself, and then drift in the energy 
will appear on this larger time scale. A meaningful long-time simulation cannot 
be expected in the unfortunate case that the timestep adaptation makes repeated 
choices in a way that resonates with some process of the system being simulated. 
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The sphere. The sphere cannot be generally uniformly subdivided into spherical 
triangles; however, a good approximately uniform grid is obtained as follows: start 
from an inscribed icosahedron which produces a uniform subdivision into twenty 
spherical isosceles triangles; these are further subdivided by halving their sides and 
joining the resulting points by short geodesies. 

Elliptic PDEs. The variational approach we have developed allows us to examine 
the multisymplectic structure of elliptic boundary value problems as well. For a 
given Lagrangian, we form the associated action function, and by computing its first 
variation, we obtain the unique multisymplectic form of the elliptic operator. The 
multisymplectic form formula contains information on how symplecticity interacts 
with spatial boundaries. In the case of two spatial dimensions, X = M 2 , Y = R 3 , 
we see that equation ( |4.36| ) gives us the conservation law 



where the vector X = {u°(j 1 V,j 1 W),u l (j 1 V,j 1 W)). 

Furthermore, using our generalized Noether theory, we may define momentum- 
mappings of the elliptic operator associated with its symmetries. It turns out that for 
important problems of spatial complexity arising in, for example, pattern formation 
systems, the covariant Noether current intrinsically contains the constrained toral 
variational principles whose solutions are the complex patterns (see Marsden and 



There is an interesting connection between our variational construction of multi- 
symplectic-momentum integrators and the finite element method (FEM) for elliptic 
boundary value problems. FEM is also a variationally derived numerical scheme, 
fundamentally differing from our approach in the following way: whereas we form 
a discrete action sum and compute its first variation to obtain the discrete Euler- 
Lagrange equations, in FEM, it is the original continuum action function which is 
used together with a projection of the fields and their variations onto appropriately 
chosen finite-dimensional spaces. One varies the projected fields and integrates such 
variations over the spatial domain to recover the discrete equations. In general, the 
two discretization schemes do not agree, but for certain classes of finite element bases 
with particular integral approximations, the resulting discrete equations match the 
discrete Euler-Lagrange equations obtained by our method, and are hence naturally 
multisymplectic. 

To illustrate this concept, we consider the Gregory and Lin method of solving 
two-point boundary value problems in optimal control. In this scheme, the discrete 
equations are obtained using a finite element method with a basis of linear inter- 
polants. Over each one-dimensional element, let N\ and N 2 be the two linear interpo- 
lating functions. As usual, we define the action function by S(q) = / Q T L(q(t),q(t))dt. 
Discretizing the interval [0, T] into N+l uniform elements, we may write the action 
with fields projected onto the linear basis as 



dWX = 



Shkoller [1997]). 




N-l 



k+l 



L({N 1( p k + N 2 (j) k+1 }, {N 1( j) k + N 2 (f> k+1 })dt. 



48 



Since the Euler-Lagrange equations are obtained by linearizing the action and hence 
the Lagrangian, and as the functions iVj are linear, one may easily check that by 
evaluating the integrals in the linearized equations using a trapezoidal rule, the 
discrete Euler-Lagrange equations given in (|3.3| ) are obtained. Thus, the Gregory 
and Lin method is actually a multisymplectic-momentum algorithm. 

Applicability to fluid problems. Fluid problems are not literally covered by the 
theory presented here because their symmetry groups (particle relabeling symme- 
tries) are not vertical. A generalization is needed to cover this case and we propose 
to work out such a generalization in a future paper, along with numerical implemen- 
tation, especially for geophysical fluid problems in which conservation laws such as 
conservation of enstrophy and Kelvin theorems more generally are quite important. 

Other types of integrators. It remains to link the approaches here with other 
types of integrators, such as volume preserving integrators (see, eg, Kang and Shang 
[1995], Quispel [1995]) and reversible integrators (see, eg, Stoffer [1995]). In partic- 
ular since volume manifolds may be regarded as multisymplectic manifolds, it seems 
reasonable that there is an interesting link. 

Constraints. One of the very nice things about the Veselov construction is the 
way it handles constraints, both theoretically and numerically (see Wendlandt and 
Marsden [1997]). For field theories one would like to have a similar theory. For 
example, it is interesting that for fluids, the incompressibility constraint can be 
expressed as a pointwise constraint on the first jet of the particle placement field, 
namely that its Jacobian be unity. When viewed this way, it appears as a holonomic 
constraint and it should be amenable to the present approach. Under reduction by 
the particle relabeling group, such a constraint of course becomes the divergence 
free constraint and one would like to understand how these constraints behave under 
both reduction and discretization. 
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